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Abstract. The adiabatic shock produced by a com- 
pact object moving supersonically relative to a gas with 
uniform entropy and no vorticity is a source of en- 
tropy gradients and vorticity. We investigate these an- 
alytically. The non-axisymmetric Rayleigh-Taylor and 
axisymmetric Kelvin-Hclmholtz linear instabilities are po- 
tential sources of destabilization of the subsonic accretion 
flow after the shock. A local Lagrangian approach is used 
in order to evaluate the efficiency of these linear insta- 
bilities. However, the conditions required for such a WKB 
type approximation are fulfilled only marginally: a quanti- 
tative estimate of their local growth rate integrated along 
a flow line shows that their growth time is at best com- 
parable to the time needed for advection onto the accre- 
tor, even at high Mach number and for a small accretor 
size. Despite this apparently low efficiency, several features 
of these mechanisms qualitatively match those observed 
in numerical simulations: in a gas with uniform entropy, 
the instability occurs only for supersonic accretors. It is 
nonaxisymmetric, and begins close to the accretor in the 
equatorial region perpendicular to the symmetry axis. The 
mechanism is more efficient for a small, highly supersonic 
accretor, and also if the shock is detached. 

We also show by a 3-D numerical simulation an ex- 
ample of unstable accretion of a subsonic flow with non- 
uniform entropy at infinity. This instability is qualitatively 
similar to the one observed in 3-D simulations of the 
Bondi-Hoyle-Lyttleton flow, although it involves neither 
a bow shock nor an accretion line. 

Key words: Accretion, accretion disks - Hydrodynamics 
- Instabilities - Shock waves - Binaries: close - X-rays: 
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1. Introduction 

The instability of the supersonic axisymmetric Bondi- 
Hoyle-Lyttleton (hereafter BHL) accretion was first dis- 
covered in 2-D numerical simulations by Matsuda et al. 

(1987) for axisymmetric accretion and by Fryxell & Taam 

(1988) , Taam & Fryxell (1989) for flows including density 
or velocity gradients. The shock surface oscillates from one 
side of the accretor to the other (so called "flip-flop" in- 
stability), leading to high-amplitude, quasi-periodic vari- 
ations of the mass accretion rate. This phenomenon was 
later confirmed by Matsuda et al. (1989, 1991, 1992), who 
showed that this process is more violent for small accre- 
tor sizes, non absorbing accretors, and high Mach num- 
bers (see also Benensohn et al. 1997, Shima et al. 1998). 
3-D numerical simulations were performed by Ishii et al. 
(1993), Ruffert (1996 and previous works for homogeneous 
media and 1997 for flows including gradients), showing 
again quasi-periodic variations of the mass accretion rate, 
although with a smaller amplitude (up to 30 per cent), and 
with deformations of the shock surface only in the imme- 
diate vicinity of the accretor. 

Livio (1992) proposed a series of possible observational 
implications of the instability. In particular, it ought to 
occur in the accretion process of a neutron star orbiting in 
a dense wind in high mass X-ray binaries (HMXB) (Taam 
et al. 1988, De Kool & Anzer 1993). It was also applied to 
the supermassive black hole SgrA* at the galactic center 
(Ruffert and Melia 1994), and even to individual galaxies 
in the intracluster gas (Balsara et al. 1994). 

Accretion onto a point like Newtonian accretor mov- 
ing supersonically in a uniform adiabatic gas is of course 
a highly idealized problem. It presents the advantage of 
depending only on two dimensionless parameters, namely 
the adiabatic index 7 of the gas and its mach number Moo 
at infinity. Although this academic formulation seems sim- 
ple, it gives rise to an instability for which no clear cri- 
terion is available yet. The extreme simplicity of this for- 
mulation leads us to expect simple laws describing the 
onset of instability, and in particular the distribution of 
timescales characterizing the instability. Numerical simu- 
lations have to include a third dimensionless parameter, 
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namely the size of the accretor in units of the accre- 
tion radius (ta = 2GM/v^ a ). Due to prohibitive com- 
putationnal cost, the smallest accretor size considered in 
3-D was r*/rA = 0.02 (Ruffert 1996 and previous works). 
We would like to be able to extrapolate the results ob- 
tained with numerical simulations to smaller accretors, 
like a weakly magnetized neutron star or a black hole, 
moving at — 1000 km s" 1 , for which r*/rA ~ 10 -5 . 

Despite numerous numerical simulations, our under- 
standing of the instability is still unsatisfactory. The in- 
stability of the accretion column for cold flows is well es- 
tablished (Cowie 1977, Soker 1990, 1991), but does not 
directly apply to the case of hydrodynamic BHL flows. A 
tentative explanation for the flip-flop instability was pro- 
posed by Livio et al. (1991), who showed that a conical 
shock becomes unstable when its opening angle exceeds a 
critical value. However, recent numerical simulations (Ruf- 
fert 1994b, 1995) suggest that the origin of the instability 
in 3-D is to be found within the subsonic flow near the 
accretor rather than at the shock surface. This lack of a 
definite physical explanation for the instability left open 
the question of whether this instability is influenced by nu- 
merical artifacts or is a natural consequence of the physics 
involved. 

A detailed understanding of the instability mechanism 
should enable us to predict how the instability is influ- 
enced by the effects of the accretor size, density and ve- 
locity gradients in the upstream flow (Ruffert & Anzer 
1995, Ruffert 1997), radiative cooling and heating of the 
gas (Blondin et al. 1990, Taam et al. 1991) and relativis- 
tic effects near the accretor (Petrich et al. 1989, Font & 
Ibanez 1998a, b), and be more confident in invoking this 
instability for the variability of accreting systems in astro- 
physics. 

To suggest such a physical mechanism is the purpose 
of the present paper. 

The local approach that we use is shortly reviewed in 
Sect. [|. A lower bound for the entropy gradient produced 
by a shock is computed in Sect. |[ We use a simplified 
formulation of the Rayleigh-Taylor (Sect. ||) and Kelvin- 
Helmholtz (Sect. ||) linear instabilities in order to estimate 
their influence on the stationary BHL flow described in 
Foglizzo & Ruffert (1997, Paper I). They are compared 
and discussed in Sect. ||. The results of new subsonic 3-D 
simulations are interpreted in the light of this analysis in 
Sect. 

2. Local stability analysis 

2.1. Local linear growth rate integrated along a flow line 

According to the numerical simulations (e.g. Ruffert 
1995), the instability of the BHL flow occurs only when a 
shock is present. The shock is a source of entropy inhomo- 
geneities and vorticity, and therefore potentially a source 
for two well known local instabilities: entropy gradients 



in a gravitational field may lead to the Rayleigh-Taylor 
instability (hereafter RTi ), and vorticity can induce the 
Kelvin-Helmholtz instability (hereafter KHi ). Note that 
entropy gradients V5 and vorticity w in the BHL flow 
are closely related by Eq. (10) in Paper I: 

w x v = TVS. (1) 

Unlike Garlick (1979) and Petterson et al. (1980) who 
used a global analysis to show the stability of the spher- 
ically symmetric accretion flows, we adopt here a local 
approach to estimate the effect of the RTi and KHi on 
the axisymmetric BHL flow. Although a global perturba- 
tive analysis would in principle lead to conclusive state- 
ments about the stability of the flow, it seems to be ex- 
cessively difficult for axisymmetric flows, where a bound- 
ary value problem must be solved in two dimensions (ra- 
dial and azimuthal) for perturbations of an incompletely 
known stationary flow (Paper I). The local approach has 
the double advantage of being mathematically tractable 
and physically understandable, although it might require 
some strong approximations. Using the same notation as 
in Paper I, flow lines are indexed by their distance w to 
the symmetry axis at infinity. We evaluate the typical lo- 
cal growth rate of each instability and integrate it over the 
time available for amplification, i.e. along a flow line w be- 
tween the shock r s h(ro) and the accretor surface r*(vj). We 
would like to check whether such a mechanism can amplify 
perturbations up to non-linear amplitudes before they are 
advected onto the surface of the accretor. We consider the 
linear growth rate |<7i(r)| of the instability (i = RT or 
KH) as obtained from a normal mode approach in an in- 
finite medium of same entropy gradient and vorticity as 
at the position r in the BHL flow. As stressed by Garlick 
(1979) in the case of spherical accretion, such a local ap- 
proximation is justified only if the distance over which the 
growth rate varies (9 log ai/dr) -1 is longer than the dis- 
tance (v/o~i) traveled during a growth time. The variation 
of the growth rate is due to the convergence and acceler- 
ation of the flow, typically on a scale r. The criterion can 
be stated quantitatively as follows: 



o~i(r)r 




d log o~i 


> 


dlogr 


V 



We estimate the quantity Ai{w) defined as: 
ft* rr*(™) \fj.( r )\ 

AM e Mr)|dt= / l -^f-dl, (3) 

Jt eh Jr Bb (-n) v { r ) 

Ai = Ma,x{Ai(zu),zu > 0}, (4) 

where we defined the elementary displacement as dl = 
vdt. Ai is a dimensionless number which depends on the 
three dimensionless parameters of the problem, namely 
the Mach number A^oo > of the flow at infinity, the 
adiabatic index 5/3 > 7 > 1 and the accretor radius r* > 
0. We aim at determining the maximum value of Ai when 
these parameters are varied. 
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The WKB approximation gives accurate results where 
its criterion (Q) is satisfied, i.e. if Ai ^> 1. A threshold 
A* exists above which the results are still significant, e.g. 
within 10%. As shown by Bender & Orszag (1978) in many 
illustrative examples of the WKB method, the WKB ap- 
proximation often gives reliable results even when its cri- 
terion is marginally satisfied. Although this may lead us to 
expect A* ~ 1, the threshold A* does not always strictly 
equal unity and naturally depends on the problem con- 
sidered. The exact determination of A* lies beyond the 
scope of this paper, and we shall assume it is of the order 
of unity. The three following situations might be encoun- 
tered: 

(i) if Ai 3> ,4*, we conclude that the criterion (0) is 
fulfilled, and that a strong linear instability is identified 
for the corresponding set of parameters. 

(ii) If Ai > A* is finite, we would conclude that a 
marginal linear instability is present, which may lead to 
a non-linear instability or saturation if the typical am- 
plitude of the perturbations is larger than exp(— Ai). For 
example, A\ ~ 3 would be enough to amplify to non-linear 
amplitudes initial perturbations of order 5%. 

(iii) If Ai < A* , the criterion (|^) is not fulfilled and our 
method does not allow us to reach a conclusive statement. 

Note that none of the 3-D numerical simulations sug- 
gest a particularly violent linear instability, since several 
accretion times ta/coo are usually needed before the insta- 
bility becomes visible (e.g. Ruffert 1995). Case (i) is there- 
fore a priori excluded in the range of parameters covered 
by these simulations, i.e. A^oo < 10, r^/r^ < 0.02. 

2. 2. Local expansion in the vicinity of a point like 
accretor 

Our local approach makes the important "continuity" as- 
sumption that the BHL flow on a Newtonian accretor of 
size r+ resembles the BHL flow on a point like Newto- 
nian accretor when r* — > 0. This allows us to make series 
expansions in the vicinity of the singularity r = 0, and 
check whether Ai diverges when r+ — > 0. Although the 
limit r — > with Newtonian gravity is unrealistic (see e.g. 
Petrich et al. 1989, Font & Ibanez 1998a, b for relativistic 
effects), it is useful in order to understand formally clas- 
sical flows, before more sophisticated effects are added. 




Fig. 1. Coefficient r\ < 1 entering the expression (g) of the 
entropy gradient created by the shock. The value of the 
adiabatic index 7 is indicated on each curve. 



the Mach number M\ associated to the velocity compo- 
nent v\±_ ahead of and perpendicular to the shock: 



2+(7-l)Ajj 
h + l)Mj 



2 1 M\ -7 + I 
7 + 1 



(5) 



Let «2_l t> e the velocity component perpendicular to and 
immediately after the shock. We write the entropy gradi- 
ent immediately after the shock as a function of Af 1 and 
its variations with respect to the curvilinear abscisse L 
along the shock, using Eq. (g): 



VS 



V 



2n v 2 c? log All 



1 V 2 ± 
27(7 



dL ' 



(2+( 7 -l)Al?)(27A4f-(7-l)) 



< 1. 



(6) 
(7) 



77 always converges to unity for large Mach numbers, when 
the kinetic energy of the gas exceeds its internal energy 
(i.e. M\ ^S> 2/(7 — 1)). The convergence is thus much 
slower for 7 close to 1 (Fig. [l]). 

According to numerical simulations, the shock distance 
seems to vary strongly from about 0.2 accretion radii for 
7 = 5/3 to apparently zero for 7 = 1.01 (r* = 0.02 in Ruf- 
fert 1996). It is not clear yet whether the shock would be 
detached for smaller accretors. Wolfson (1977) remarked 
that for 7 close to one, energy is soaked up by the internal 
degrees of freedom of the gas, therefore not contributing 
to support the shock through the kinetic pressure, thus 
favouring an attached shock. This leads us to consider 
successively the cases of attached and detached shock. 



3. Entropy distribution produced by the shock 

3.1. Entropy gradient along the shock 

For the sake of simplicity, in what follows we choose the 
same units as in Paper I, such that the ratio of the mean 
molecular weight to the gas constant n/lZ = 1, without 
loosing any generality. The Rankine-Hugoniot conditions 
(e.g. Landau & Lifshitz 1987) imply that the entropy jump 
AS across an adiabatic shock is an increasing function of 



3.2. Case of a shock attached to a point like accretor 

We deduce from Eqs. (El), (E2), (E4) in Paper I that 
along a shock attached to a point like accretor with an 
angle s h, for r <C r\, the velocity scales like: 



via 



vi\\ 



/2GMV S 



L J 
( 2GM \ ' 



cos ■ 



V L 



2 



7 + 1 



V2\\- 



V2A 



(8) 
(9) 
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Using Eq. (0), the entropy gradient immediately after the 
shock, close to the point like accretor is: 



|V5| 



7-1 



1 



7+1 
7- I 



tan 



Csh 

T 



i 

T 



(10) 



The entropy produced by the shock is a decreasing func- 
tion of L if the shock is attached. It depends on the Mach 
number only through the shock opening angle # s h- 

3. 3. Case of a detached axisymmetric shock 

Let r s h(0) be the shape of the axisymmetric shock surface, 
using polar coordinates (r, 9) centred on the accretor. Let 
rv = r s h(7r) be the distance of the shock from the accre- 
tor, along the symmetry axis. Eq. (|^) indicates that when 
M-oo 3> 1) the entropy jump along the symmetry axis is: 

ASK) ~ — ^-log-Mi > -H-logA^ . (11) 
7 — 1 7 — 1 

Far from the accretor, the shock surface approaches the 
Mach cone of semi-angle 9 S defined by sin# s = I /Moo 
{Mi ~ 1). The entropy jump therefore decreases from 
A5(tv) ahead of the accretor to zero far from it. Since 
the entropy gradient immediately after the shock surface 
vanishes both on the symmetry axis and far from the ac- 
cretor, the maximum entropy gradient VS| max is reached 
on a circle r s h(5 max ) corresponding to an intermediate az- 
imuthal angle max . 

Let us denote by a the angle between the flow line and 
the vector perpendicular to the shock surface, before the 
shock, so that i>u_ = v\ cos a. Defining the dimensionless 
function £ along the shock surface r s i 1 {9) as 

rv 91ogA^i 



C(r*(0)) 
Eq. 



cos a dL 
is rewritten as follows 



<Cn 



2 V 



1 < 



7" 

V2± 



i'T 



- 1 sin 2 



(j+l)Ml 7+1 



< 



2 + { 1 -l)M\ 7-1 



(12) 

(13) 
(14) 



The Rankine-Hugoniot condition (Eq. [li]) and Eq. ( |l3| ) 
provide us with both a lower and an upper bound for the 
maximum entropy gradient produced by a strong shock 
{r, ~ 1 for Moo » 1): 

2(7+l)Cmax ^ lTTOI x 2 ( n 



(7 -I) 2 



> IV5I. 



> 



7' 



1 



(15) 



The value of the maximum Cmax of the function £ depends 
only on the properties of the supersonic flow before the 
shock, and on the shape of the shock surface r s h(6>). The 
function C*(?+) is defined as the minimum value of Cmax 
for all possible continuous mathematical curves r s h(0) sat- 
isfying r sh (7r) = ?v and dr sh /d9(ir) = 0: 

C*(fv) = Min {Cmax, any curve r sh (6) , r sh (7r) = ?v} . (16) 
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Fig. 2. The coordinates of the accretor are (0,0) on the 
upper plot. The hyperbolic trajectories of the gas ema- 
nating from x < are represented by dotted lines. Four 
particular shock shapes are plotted for r v =0.1 accretion 
radii. The short dash line is orthogonal to the symmmetry 
axis, and the long dash curve is orthogonal to the super- 
sonic flow lines. The curves with circles are the optimal 
shapes leading to the smallest entropy gradients corre- 
sponding to polynomial shapes of the variable y 2 (empty 
circles) and y 1 ^ 2 (filled circles). The value of £ along these 
four curves is displayed on the bottom plot. 

Thus we obtain a lower bound for the maximum entropy 
gradient iVSJmax produced by a detached shock standing 
at the distance rv from the accretor: 



VSI > 

I max — 



2 C* 



7 



1 TV 



(17) 



In the framework of the approximation of the supersonic 
trajectories by hyperbolae, for Moo 3> 1, C* depends only 
on the distance r n of the shock. We have computed nu- 
merically the function £*(rv) using a polynomial approxi- 
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Fig. 3. Coefficient characterizing the minimum value 
of the maximum entropy gradient immediately after an 
axisymmetric shock, depending on its distance rv to the 
accretor. The shock shapes are the same as in Fig. ||. 



mation of the shock shape and a downhill simplex method 
for the minimization (Press et al. 1992). Powell's method 
was also used with comparable results. Satisfactory re- 
sults were obtained with a polynomial x(y) of order 3 in 
y 2 . The overall minimum seems to be reached by the sin- 
gular curve scaling like x + r„ ~ ?/ 3 / 2 near the symme- 
try axis, which can be approached by a series of regu- 
lar polynomials (Figs || and ||). For comparison, a plane 
shock orthogonal to the symmetry axis produces typically 
twice as much entropy gradients than the minimum value 
(Cmax/C* ~ 2). More realistic is the curve orthogonal to 
the supersonic flow lines (a = 0) , which produces entropy 
gradients about 20% stronger than the absolute minimal 
value. According to Fig. ||, the coefficient (* > 0.1 for real- 
istic shock distances, i.e. r T < r&- The maximum entropy 
gradient along the shock stands near L ~ r n (Fig. ||). 

Fig. H shows that the mathematical curves approach- 
ing the minimum value of the entropy gradient after the 
shock are not very different from the physical shock shapes 
observed in numerical simulations (e.g. Fig. 0). The lower 
bound C* might therefore be a good approximation of the 
realistic value of Cmax, within a factor two. 



3-4- Entropy gradient in the subsonic flow between a 
detached shock and the accretor 

If the boundary condition on the surface of the accretor 
allows a high enough mass accretion rate, the maximum 
entropy gradient immediately after the shock corresponds 
to a flow line tu max converging to the accretor. With the 
entropy remaining constant along each flow line after the 
shock, the gradient of entropy across the flow is simply the 
gradient immediately after the shock surface modified by 
a geometrical factor. If the convergence towards the ac- 
cretor were along straight flow lines, the entropy gradient 
would simply increase like 1/r. However, flow lines are not 
radial, and we can take this into account by including a 



geometrical factor 6(r), so that the distance between two 
neighbouring flow lines converging towards the accretor 
scales like r/S(r): 



5{T) . 

r s h (VAjsh 



(18) 



where r(0) is the shape of the flow line w max , and (VS%h 
is the entropy gradient on this flow line immediately after 
the shock, at a distance r s h, so that <5(r s h) = 1. We made 
a distinction in Sect. 4 of Paper I between the directions of 
regular and singular accretion on a point like accretor. The 
distance between adjacent accreted flow lines decreases 
like ~ r when r — > in a direction of regular accretion 
(lim r ^o S(r) is finite), whereas it decreases much faster 
in a direction of singular accretion (lim r ^o S(r) diverges). 
Since we proved in Paper I that accretion for a gas with 
7 = 5/3 is always regular, we know that the instability of 
the BHL flow does not rely on the presence of directions of 
singular accretion. We can therefore assume that lim r ^o 8 
is finite in our analysis of the instability. 

We obtain from Eqs. ( ]lq ) and ([18j) upper and lower 
bounds for the entropy gradient between the detached 
shock and the accretor, along the flow line n7 max , at high 
Mach number: 



2 (7 + l)Cmax r sh S 2Cmax ^sh $ 

( 7 - 1)2 r T r - 1 !max - 7 „ 1 ^ r 



(19) 



Note that the ratio r s h/?V is of the order of 2 2 in Fig. ||. 

4. Rayleigh Taylor instability 

4-.1. A simplified formulation of the RT instability 

Let us consider a stratified gas in a gravitation field. An 
entropy gradient can act either in a stabilizing or desta- 
bilizing manner depending on whether it does or does not 
contribute to support the gas against gravity. The verti- 
cal gravity field g and pressure forces at equilibrium are 
related as follows: 



1 

Po 



-VP 



(20) 



The Brunt-Vaisala, frequency N is the frequency of oscilla- 
tions of the stratified gas, in the limit of small wavele ngth s 
perpendicular to the gravity field (see Appendix Al). 
<7rt = (—A 2 ) 2 is the local growth rate of the RTi when 
the entropy decreases upwards: 



o&r = -N 2 = 1 — -g ■ VS* . 

7 



(21) 



If the flow is sheared with height, perturbations with a 
short horizontal wavelength in the direction perpendicu- 
lar to t he fl ow velocity grow at the same rate (see Ap- 
pendix A..1). 
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Fig. 4. Effective gravity, in GM/r 2 units, in the case of 
spherical accretion for various values of 7. 



4-2. Effective gravity in the comoving frame 

From the principle of equivalence, there would be no RTi 
in a gas falling freely in a gravitational potential, because 
the effective gravity would then be zero. So in a reference 
frame falling with the gas, the effective gravity g c g driving 
the RTi is opposite to the pressure force: 



ffeff 



p 



(22) 



According to Sect. 4.6 in Paper I, the pressure close to 
the accretor is spherically symmetric to first order for 
7 = 5/3. This leads us to neglect the negative contri- 
bution of the azimuthal pressure force to the scalar prod- 
uct VP • VS in Eq. The radial pressure support 
decreases from the subsonic region to the supersonic re- 
gion. The effective gravity, calculated analytically in the 
case of spherical accretion, is displayed in Fig. ^ for vari- 
ous adiabatic indices. In the subsonic region, the effective 
gravity is more than 50% of the gravity of the accretor 
for any value 1 < 7 < 5/3. The best pressure support 
is reached, of course, for 7 = 5/3, for which the effec- 
tive gravity in the subsonic region is at least 75% of the 
gravity. We shall assume that the effective gravity for an 
axisymmetric flow is comparable to the effective gravity 
in the spherical case, thus constraining the dimensionless 
gravity parameter g* G [0.5, 1] in the subsonic region of 
the flow: 



ffcff 



GM 



(23) 



Jj.,3. RTi efficiency in the BHL flow 



4.3.1. General expression of the RT efficiency 

According to Eq. (^l|) , the local growth rate of nonaxisym- 
metric perturbations with a short wavelength perpendicu- 
lar to both the flow lines and the radial direction (similar 
to the case k x — ► 00, k v = in Appendix A.l) is approxi- 
mated as follows: 



2 7-I * GM . 
^rt = 9 — 2- V5 sin/3 

7 



(24) 



where (3 is the angle between the flow line and the radial 
direction: 



tan p = — . 



(25) 



The perpendicular wavelength of a non axisymmetric per- 
turbation decreases geometrically as the flow is advected, 
so that the growth rate of the RTi stays maximum. 
Defining the free fall velocity vg as 



^2GM ^ 2 



(26) 



the integrated efficiency of the RTi defined by Eq. 
follows from Eqs. © and ©: 



7-1 
2 7 

r M ro ) lV(r iid/ 
ff *^( r |VS|) 5 sinM- 

r sh M ' v r 



(27) 



With the estimates of t he e ntropy gradient and effective 
gravity of Sects. || and we are now able to evaluate 
Eq. ( p7| ) for the different topologies discussed in Paper I. 
Since the entropy along the shock decreases away from 
the symmetry axis, the stratification is potentially linearly 
unstable only in the region where (3 > 0. The conservation 
of angular momentum implies that (3 > in the supersonic 
flow before the shock. Since (3 is likely to increase across 
the shock, we conclude that the stratification is locally 
unstable immediately after the shock surface. We denote 
by rp (9) > the surface where the velocity is radial, thus 
delimiting a region of unstable stratification. 

4.3.2. Attached shock 

Using Eqs. (|l(]), ( p4| ) close to accretor (?y ~ 1), and the 
relation sin/3 = «2_l/^2, the local growth rate is: 



g \ 2 v s 



(28) 



This growth rate must be integrated along a flow line w 
between the azimuthal angles 8 s h and 9 so corresponding to 
the shock and the sonic surfaces respectively. The length 
of this path of integration is of the order of r(6 s h — 9 so ). 
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Estimating the velocity after the shock from Eqs. (||) and 
(§) gives: 



•Art (a?) 



CRT 



dl. 



7+1 
2 7 



7sh 



sin 



2 0«1 



(7-1) 2 
47 



(29) 

- < 23.(30) 



The efficiency depends strongly on the azimuthal size of 
the subsonic region reaching the accretor (6* s h — SO )- This 
parameter is unfortunately unknown, and we only obtain 
an upper bound on the efficiency of the RTi : the timescale 
of the RTi instability is at best comparable to the advec- 
tion timescale for a shock attached to the accretor, how- 
ever small the accretor might be. 

4.3.3. Region of supersonic accretion near a point like 
accretor with a detached shock, 7 < 5/3 

The sign of (3 in the supersonic region might simply pre- 
clude the instability (/3 < for 7 = 5/3). Let us show 
that even if the flow lines were bent in the unstable direc- 
tion [p > for 7 ~ 1), the RTi would become negligible 
when the gas approaches a point like accretor. For this we 
wish to check that the divergence of the entropy gradient 
when r — > is not fast enough to make the RTi growth 
time shorter than the free fall time. Using Eq. ([ljj]) and 
Eq. (|J) , the radial dependence of the RTi growth rate in 
the supersonic region scales like: 



(31) 
(32) 



where we have used Eq. (56) of Paper I for an upper bound 
of/3(r — > 0) (d log (3/d log r > (5— 3 7 )/2), and the decrease 
of the effective gravity. Using the free fall approximation 
of the velocity, the contribution of the region surrounding 
a point like accretor to the efficiency of the RTi scales like: 




CRT 



dr 



cos/3 



/ 5-3-1 \ 



(33) 



From the convergence of the integral when r — - » we 
deduce that the local growth time is much longer than 
the free fall time, and the RTi can be neglected there if 
7 < 5/3. 

4.3.4. Region of subsonic accretion along the maximum 
entropy gradient 

Since the RTi is driven by the entropy gradient, it is nat- 
ural to evaluate a lower bound for its efficiency along the 
flow line w max associated with the maximum entropy gra- 
dient, at high Mach number (77 ~ 1), using Eqs. (|19| ) and 



Art > Min^RT, 

Cm ax fsh. 



(34) 



Min _4 RT = 



7 



1 dl 



— (5g* sin ft)?— . (35) 



The Bernoulli equation can be written in terms of the 
free fall velocity, and approximated for M 2 ^ 3> 2/(7 — 1), 
inside the sonic radius: 



1 



(7 - l)M 2 



= 4 



1 



(Tf-l)Ml 



,(36) 
(37) 



where we have neglected Vqc compared to vg between the 
shock and the accretor, since the shock distance is typi- 
cally shorter than the accretion radius («gf(rA) = «oo)- The 
ratio Vff/v decreases in the subsonic region between the 
shock (r s h) and the sonic point (r son i c ). Applying Eq. fl37| ) 
at the sonic point (M = 1) and at the shock (M ~ M s h), 
we obtain the following range: 



7' 



1 



7 



1 



Ug 7 + 1 



7- 



1 



(38) 



If 7 < 5/3, the flow line tn max reaches the accretor in 
the supersonic region, and we know from the preceding 
section that the contribution of this region to the integral 
is negligible. 

Let us examine each of the terms of Eq. (^) for 
7 = 5/3, which is supposed to be the most unstable case 
according to numerical simulations. 

(i) 0.75 < g* < 1 according to Fig. [|, and thus we 
estimate g* ~ 0.9, 

(ii) the geometrical factor 6 is finite since accretion 
with 7 = 5/3 is always regular (Paper I), and is assumed 
to be of the order of unity, 

(hi) 2 < vs/v < 4 according to Eq. (|3§|), and thus we 
estimate Vg/v ~ 3. 

Because these three contributions to the integral are 
finite, we replace each of them by their mean value and 



approximate Eq. (35) as follows 



Min ^ RT ~ 1.0 



Cn 



0.15 2*rv 



" . x a dl 
sin 2 p — 

T 



(39) 



We showed in Paper I that (3 < along the sonic surface 
for 7 = 5/3, thus rp{vu) > can be estimated as a frac- 
tion of the shock distance r w . This precludes the possibility 
that the integral in Eq. (^9|) might diverge when the size 
of the accretor is decreased. This further suggest that the 
efficiency of the RTi instability should not increase much 
when the accretor size is much smaller than the shock dis- 
tance (r* <C r w ). In order to estimate this efficiency, let 
us remark from Fig. [?] that sin (3 is close to unity immedi- 
ately after the shock (see also numerical simulations as in 
Fig.^). Since we expect that (3 > over a sizable fraction 
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Fig. 5. Typical dependence of the Kelvin-Helmholtz insta- 
bility on the wavelength A of the perturbation, described 
by the function \ (full line). The dotted line shows the 
slope at long wavelengths. 

of the shock distance, this integral is likely to be of or- 
der unity. According to the typical values of £max deduced 
from Fig. ||, we conclude from Eq. (|3|) that Min „4 RT is of 
the order of unity, hardly more. Since the lower and upper 
bounds in Eq. differ by a factor (7 + l)/(7 — 1), and 
since the RTi growth rate scales like the square root of the 
entropy gradient, we finally estimate for 7 = 5/3: 



1.0 < „4 RT ^ 2 -° 



(40) 



Thus the integrated efficiency of the linear RTi does not 
diverge, even in the case of a point like accretor moving 
at high Mach number in a gas with 7 = 5/3. Although 
small, this efficiency of order unity is not negligible. 

5. Kelvin Helmholtz instability 

5.1. A simplified formulation of the KH instability 

The maximum linear growth rate of the KHi in an inviscid 
fluid is comparable to the maximal vorticity of the flow 
in equilibrium. From the studies of the instability with 
various velocity profiles (see an overview in Drazin and 
Reid, 1981), one can extrapolate the following growth rate 
ctkh and optimal wavelength A max of the instability cor- 
responding to a vorticity profile with a maximum |w max | 
and a gradient width h: 



|ckh| = a |w max | x(-y^), with a « 0.2, 



with 



(41) 
(42) 



where x( x ) x(l) = 1 is a function of the wavelength A of 
the perturbation. Its typical shape, obtained by a numeri- 
cal solution of the Orr-Sommerfeld equation, is plotted in 
Fig. |. 

The effect of compressibility on the linear instability 
can be studied by solving the linearized equations for a 
sheared plane flow (see Appendix |A.2j ) , with the following 
velocity profile: 

v(z) = ^tanh-, (43) 



tanh — 
2 h 




Fig. 6. Maximum growth rate of the linear KHi and as- 
sociated wavelength, as a function of the Mach number 
A^kh = vq/cq for a uniform sound velocity cq, and a ve- 
locity profile given by Eq. (|43|). The growth rate ctkh is 
displayed in units of the vorticity maximum w max = v /h 
(solid line). The dimensionless parameter k measures the 
optimal wavelength A max in units of the width of the vor- 
ticity peak h (dashed line). The rigid boundaries are at 
±5/i. Within the accuracy of our numerical relaxation 
method, these curves are independent of the adiabatic in- 
dex in the range of interest 1 < 7 < 5/3. 

Neither the value of 7, within the range [1,5/3], nor the 
presence of an entropy gradient across the flow influences 
the growth rate and the wavelength of the most unstable 
mode by more than a few percent. 

A determining quantity for the KHi in compressible 
flows is the relative Mach number A^kh measured in the 
frame comoving with the flow line of maximum vorticity: 



_ \v(h/2) -v(-h/2)\ hw n 



CO 



(45) 



(44) 



with Co being the sound speed on the line of maximum 
vorticity (z — 0), and v(z) the z-dependent velocity of the 
flow in the y-direction. The growth rate decreases by a 
factor of 2 between the very subsonic flow (which mimics 
the incompressible case) and A^kh = 1-3, as shown on 
Fig. [| So the difference of velocities within the vorticity 
peak must be less than the sound speed (i.e. Mkh < 1) in 
order to allow pressure forces to act in the KHi mechanism 
as efficiently as in the incompressible limit. Note that the 
limiting values of ctkh and k for A^kh — ► in Fig. ^| agree 
with the values of a and k stated in Eqs. (|4l]) and (f42|). 

5.2. Application to the BHL flow with a detached shock 

According to Eq. ([!]) , the same symmetry arguments used 
for the entropy gradients in Sect. [| apply to the vorticity. 
A surface of steepest increase of the vorticity therefore 
connects the shock to the accretor, where the vorticity 
diverges for r — > (Eqs. 49-59 in Paper I). If this surface of 
maximum vorticity is not too different from a flow surface, 
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Fig. 7. Lines of constant entropy superimposed on a 
greyscale map of the vorticity in a numerical simulation 
of BHL with 7 = 5/3 and Moo = 10 (model FS in Ruf- 
fert 1994b). In a good approximation, the line of steepest 
increase of the entropy gradient coincide with the line of 
steepest increase of the vorticity, and corresponds to a line 
of flow. 

the KHi will develop most efficiently along these particular 
flow lines. Numerical simulations, illustrated by Fig. [?], 
together with Eq. (|l|), suggest that the surface of steepest 
increase of the entropy gradient and the surface of steepest 
increase of the vorticity coincide with the flow line denoted 
by tn max , such that #o( ro max) is a direction of maximum 
entropy gradient and vorticity. 

We use Eqs. (Q) and ( |l9| ) along these flow lines to write 
the vorticity as follows: 



2Cmax <j r sh V 

7(7 — 1) M 2 r n r 



> 



(46) 
(47) 



6 is related to the width h(r) of the vorticity peak (h s h = 
M r sh)), and thus to the optimal wavelength of the KHi 
through Eqs. (ph and @: 



S(r) 



W) 



r 



(48) 



Although the vorticity defined in Eq. ( f47| ) may become 
arbitrarily large when 7 is close to unity, the time-scale 
of the KHi is limited by compressibility effects. Let us 
estimate the Mach number Mkr defined in Eq. (p)|), using 
Eqs. © and ©: 



M 



KH 



2Cn 



hsh 1 



7(7 - 1) rv M' 



(49) 



Paradoxically, Eq. (|49|) indicates that compressibility ef- 
fects are stronger in subsonic regions. We introduce in 
Eq. ( f49| ) the minimum value A4 s h of the Mach number 
after the shock, defined by the Rankine-Hugoniot jump 
conditions: 



1- 



1 



M 



2 7 
0.4 



KH 



7 

0.6 



(7- 



1) 



S,max "sh 

015 77 



M s 



M 



0.15 



M : 



(50) 
(51) 
(52) 



where Eq. ( |52| ) assumes 7 = 5/3. Noting that the Mach 
number increases from the shock to the sonic surface 
(Ai > A4 s h), and that h^/r^ is of the order of unity, 
we conclude that the effect of compressibility on the KHi 
can be neglected for 7 — 5/3. From Eq. ( f4l| ) and (|47|), we 
can write the minimum KHi growth rate, at high mach 
number (./Woo 3> 1, V ~ 1), as follows: 



okh 



> 



2«Cn 



X 



7(7 - 1) M 2 r. 



r s h v 
r 



(53) 



7T 



Let us estimate the minimum efficiency Min „4kh of the 
KHi mechanism along the flow line w m:ix : 



Min i K H(Wmax) 



r s h 



7(7 - 1) r n 







X 2 r 



(54) 



The width /i of the entropy gradient decreases geometri- 
cally with r, and the optimal wavelength of the KHi de- 
creases according to Eq. (fl2|). By contrast, the wavelength 
A, parallel to the flow line, of a Lagrangian perturbation 
must increase when advected, since the flow is accelerated. 
It is therefore not possible to maintain the KHi with its 
local maximum growth rate during the advection of the 
perturbation. 

If we denote by A; = nhi <c r s h the initial wavelength 
of the perturbation, Fig. [| indicates that it becomes un- 
stable only when advected towards a region where the 
gradient width is h(r) < 2X/k. We rewrite Eq. (|54|) using 
Eq. ©: 



Min .4 K h 





r 


\x(h/hiY 


1 ^sh d/ 


7(7 - 1) »v 




h/hi 


M 2 hi r sh 



.(55) 



Neglecting the increase of A^ due to the acceleration com- 
pared to the linear decrease of h(r), the ratio A max (r) /\ = 
h(r)/hi decreases linearly to zero when r — > 0. The inte- 
gral in Eq. (Eq) is approximated for — > by extracting 
some average values 8 and M. from it, and integrating the 
function \ described in Fig. |[ 





WW] 


Jo 


h/hi 



1 h s h dl 

M 2 hi r sh 



6 

M 2 Jo 
8 
M 2 



dx. 



2.1 



(56) 
(57) 
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Fig. 8. Splitting of the entropy gradient peak generated 
by the shock into several narrower peaks. 



Min Akh 



0.4 



( 7 -l) 2 0.2 0.15 2irv V M 

Introducing 7 = 5/3 into Eq. d58|), and approximatin 
M ~ (M s h + l)/2 we obtain: 

Min^H-O.S-^^-? 



0.2 0.15 25rv ' 

We deduce from Eqs. (|l9| ) and ( ^9| ) an estimate of the KHi 
efficiency for 7 = 5/3: 



0.5 <, Akh & 2.0 



In particular, Akh does not diverge when r* — » 0. Ac- 
cording to Eq. ([56"|), the integral in Eq. ( |55| ) is regular at 
r = 0, thus its value should not depend on the accretor 
size if r+ <C r n . Consequently, the KHi efficiency should be 
little influenced by the size of the accretor if r* Cr,. The 
above calculation indicates that a whole range of wave- 
lengths leads to comparable efficiencies if r* <C r x . Al- 
though the integrated efficiency remains limited, the vari- 
ability timescale is directly related to the wavelength of 
the perturbation. Numerical simulations show that the 
gradient width, a priori comparable to the radius of cur- 
vature of the shock, is split into several maxima of much 
smaller width, typically of 5-10 degrees in a snapshot of 
the flow before instability (Fig. ||). One may imagine that 
the gradient width decreases with time, as more and more 
gas is accumulated from the downstream side of the ac- 
cretor, until the gradient width is short enough to start 
the KHi . We must also note that the generation of vor- 
ticity at the interface between the grids in the multi-grid 
PPM numerical technique might influence the distribution 
of vorticity in the flow. 



6. Discussion of the efficiencies of the two 
instabilities 

The local time scale of the instability in the subsonic re- 
gion scales like the advection time onto the accretor, so 
that the efficiency integrated along a flow line is always of 
order unity (Eqs. ^ and ^30|). Assuming that the thresh- 
old of the WKB approximation is also of the order unity 
(A* ~ 1) and following the statements of Sect. |^, we stand 
in the uncomfortable position between case (ii) and case 
(iii). Despite this uncertainty, it is interesting to compare 
the analyze the results of numerical simulations in the 
light of these physical mechanisms. 

The efficiency of the instabilities naturally increases 
with the strength of the shock, as can be seen from the 
function 77 (.Mi) (Eq. |?] and Fig. Ilj). Numerical simulations 
with 7 = 4/3 (Ruffert 1995) and 7 = 5/3 (Ruffert & Ar- 
nett 1994, Ruffert 1994b) show the same trend: the flow 
is stable with Moo = 1-4 and unstable with Moo = 3 and 
10, with detached shocks in all these situations. The high- 
est entropy gradients are obtained for M 2 ^ > 2/(7 — 1), 
thus requiring higher Mach numbers for nearly isother- 
mal flows. Increasing the Mach number does not increase 
the efficiency of the instabilities indefinitely: our estimates 
show that this efficiency saturates when the increase of the 
entropy gradient is compensated by the decrease of the 
advection timescale. We see no obvious justification for 
a possible divergence of any of the numerous dimension- 
less parameters introduced (£, 5, g* , a) if the Mach number 
(59) Moo is increased, or if the accretor size is decreased, and 
therefore expect the integrated efficiencies to be smaller 
than 2. 

Our calculations seem to indicate that the efficiency 
should increase when 7 approaches unity (Eqs. [35, 38 and 
|58|), whereas the most unstable flows observed in simula- 
tions correspond to 7 = 4/3 and 7 = 5/3. Nevertheless, we 



(58) 



(60) 



showed in Sect. 4.3.2 that the efficiency of the RTi is finite 
for 7=1 when the shock is attached to the accretor. Thus 
this paradox would disappear if a critical adiabatic index 
exists below which the shock is attached to the accretor, 
as suggested by Wolfson (1977). 

The axisymmetric KHi must be considered together 
with the effect of stratification, which can be stabilizing 
or destabilizing depending on the sign of (3. The influence 
of stratification on the KHi can be estimated quantita- 
tively by comparing the time scales associated to each 
physical process, in the same spirit as was done with the 
Richardson number. This ratio informs us about a possi- 
ble stabilization of the KHi by buoyancy forces if (3 < 0, 
or which of the two instabilities is the fastest if (3 > 0. The 
ratio of the timescales can be deduced from Eqs. J24|), (41) 
and (Eo) as follows: 



CRT 

ckh 



7(7 - 1) 



, y2 Uff sin 2 (3 
M — (rVS n 



<r*-(6i) 
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We estimate this ratio both immediately after the shock 
and at the sonic point. At the shock, the effective gravity 
is strong (g* ~ 1), and Eqs (pit), ( pl[ ) and (p0|) provide us 
with a lower bound for Eq. (j6l|): 




(rsh) > 3.8(7 - 1) 



CRT 
0"KH 

Applying this equation to 7 = 5/3, we obtain: 
sin5 $ ( 0.15 V 0.2 ( 2%r^ 



(62) 



CRT 
CKH 



(r«h) > 2.6 



X 



6(n 



a 



(63) 

Since (3 is large immediately after the shock, Eq. (|63| ) sug- 
gests that the RTi is more unstable than the KHi there. 
At the sonic point (M = 1), using Eqs. (|l|) and @), 
Eq. ( |6l| ) becomes: 



(rsonic) > 6.671(7-1) 



sin 2 /3 



CRT 
CKH 



Applying this equation to 7 = 5/3, we obtain: 
0-rt , > - _ sin 2 " 

I "sonic) 
CKH X 




(64) 



0.75 



0.15 

^Cmax 



(L2 / 22^ 



(65) 



Let us first remark from Eq. (71) of Paper I that if 7 = 5/3, 
< lim r ^o(<91og/3/<91ogr) < 1, indicating that the con- 
vergence of sin 2 f3 towards zero is slower than that of 
x(r) ~ r/2 close to a point like accretor. Thus the KHi 
ultimately becomes stabilized by buoyancy forces near a 
point like accretor, when the wavelength of the perturba- 
tion is much longer than the vorticity gradient width. 

According to Eq. (^59), a perturbation with a wave- 
length such that x = 1 a t the sonic point would be lo- 
cally unstable to the KHi with hardly any stabilization by 
the buoyancy forces if \/3\ < 1 degree. The angle (3 at the 
sonic point is negative for 7 = 5/3 (see Paper I) and much 
smaller than at the shock, especially if the sonic surface 
is close to the accretor. We know from Paper I that the 
sonic surface is attached to the accretor if 7 = 5/3. Thus 
the flow line reaching a point like accretor along the sonic 
surface is unstable to the KHi without being stabilized 
by the buoyancy forces. The shortest timescale associated 
to the KHi therefore depends on the size of the accretor, 
while the longest is directly related to the shock distance. 

Although the RTi is dominant near the shock, and the 
KHi might prevail closer to the sonic surface, it is hard to 



disentangle the two mechanisms, which act in a combined 
way when the shock is detached. The situation is different 
when the shock is attached to the accretor, since there 
is no vorticity maximum in this case. It is interesting to 
note that the linear efficiency estimated in Eq. ( |30| ) for 
an attached shock is not much smaller than the efficien- 
cies derived for a detached shock. Numerical simulations 
indicate that accretion flows with an attached shock are 
generally more stable than those with a detached shock 
(Ruffert 1995, Zarinelli et al. 1995). In the simulations by 
Ruffert (1996) with 7 = 1.01, the instability is neverthe- 
less clearly observed, producing fluctuations of the mass 
accretion rate of about 7% (simulations with 7 = 4/3 and 
5/3 produce typical fluctuations of 14% and 23% respec- 
tively). This suggests that the RTi alone may destabilize 
linearly the flow when the shock is attached, while the 
combined action of the two mechanisms when the shock 
is detached is more efficient and leads to higher non linear 
amplitudes. 

7. Instabilities in simulations with subsonic flows 

In order to separate more clearly the different effects that 
might contribute to the formation of the instabilities, we 
performed numerical simulations of accretion from a flow 
with subsonic bulk velocities at infinity, but with a gradi- 
ent in entropy. The rationale is the following. When the 
accretor moves supersonically relative to a surrounding 
medium a bow shock will form. The shock transforms the 
initially (at infinity) homogeneous medium moving at su- 
personic velocities into a subsonic flow with an entropy 
decreasing away from the axis of symmetry. 

We did a set of simulations to mimic these conditions, 
of which we will present here the most important results 
relevant to the topic of this paper. The accretor radius is 
chosen to be = 0.02rA and the same absorbing bound- 
ary conditions are used as in Ruffert (1994a). The simu- 
lations are done on 7 nested cartesian grids, the zone size 
on each finer grid being a factor of two smaller than of the 
next coarser grid. The largest grid spans 8 accretion radii. 
Matter is evolved hydrodynamically using the "Piecewise 
Parabolic Method" , assuming it is a polytropic gas with 
an adiabatic index of 5/3. 

The flow at infinity is in direction of the positive x-axis, 
with a constant velocity along x. The pressure is uniform 
at infinity. In order to maintain pressure equilibrium, the 
maximum of entropy AS at y — at infinity along the x- 
axis is offset by a minimum in density with the following 
shape: 



Poo{x,y,z) 

= 1 — e exp 

Po 

with £ = 1 — exp 



2 

7-1 



AS 



(66) 
(67) 



The entropy difference was chosen to be AS — 2 (e — 
0.55, rAV5 max ~ 1.33) which is comparable to the value 
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Fig. 9. Contour plot of the density distribution together with the instantaneous flow velocities. The contours are 
spaced logarithmically with intervals of 0.1 dex, the contour corresponding to 10° is annotated. The time elapsed since 
the beginning of the simulation together with the velocity arrow unit is shown at the top right. The accretor is at the 
center with coordinates x = 0, y = 0. The right figure is a zoom of the left figure enlarged around the accretor. 



observed in the numerical simulations, e.g. in Fig. 17d of 
Ruffert (1994b) and Fig. fj] of this paper. The velocity at 
infinity was chosen constant, such that the Mach number 
at infinity is 0.6 far from the axis. The density profile 
(Eq. |66|) and the uniform pressure imply that the Mach 
number decreases by a factor (1 — e) 1 / 2 along the axis. 

Fig. ^| shows a snapshot of the density and velocity dis- 
tribution of such a model. At the point in time at which 
the snapshot was taken two transient vortices are promi- 
nent at (x ~ —0.1, y w ±0.15), corresponding to a buoyant 
vortex ring rising against the flow. Such rings were already 
observed in numerical simulations by Koide, Matsuda & 
Shima (1991, Fig. 12). When boundary conditions were 
chosen for which there is little or no accretion at the sur- 
face of the star, e.g. by Shima et al. (1985), Fryxell, Taam 
& McMillan (1987), Matsuda et al. (1991), such vortex 
structures appeared too. 

A movie showing the temporal evolution of the den- 
sity distribution shows that the flow is not stationary. Just 
from an inspection of the contour plots and movies no ob- 
vious difference is apparent between the instabilities gen- 
erated in these simulations and the instabilities present in 
the simulations of three-dimensional BHL-accretion (e.g. 
Ruffert 1996, and references therein). 

Since the entropy gradients are present over the whole 
domain of the simulation, it is important to check the con- 
tribution of the region far from the accretor (e.g. r > ta) 
to the efficiency of the instabilities. Thus we decompose 
the global efficiency of each instability into the contribu- 



tions Ai and A°° defined as follows: 



^dL 



(68) 
(69) 



Far ahead from the accretor, the trajectories are nearly 
parallel to the symmetry axis, and the growth rate of the 
RTi along a flow line indexed by w ~ ta can be estimated 
from Eq. (pj) with sin/3 ~ vj/r: 



CRT 



7-1 — q GM 
vj\7b 



A 



RT 



< 



7 

7-1 
27 



3 J 
T 2 



r A VS 



du 



(1+m 2 )S 



< 0.6, 



(70) 

(71) 
(72) 



where we have used the parameters of the simulation, 
namely r^/rx — 8, and 7 = 5/3. The convergence of 
the integral when — > 00 ensures that the main con- 
tribution of the RTi comes from the region close to the 
accretor, within a few accretion radii. Nevertheless, the 
value of A^ T might not be fully negligible compared to 
Art as estimated in Eq. (p0|) . 

An important difference between this model and the 
supersonic models with a bow shock is that here the Ber- 
noulli constant B{r, 9), defined by B = v^, /2 + c 2 ^ / (7 - 1) 
is not uniform far away from the accretor, since the sound 
speed is not uniform at infinity. The non-uniformity of B 
adds an additionnal term in Eq. (|l|) (Eq. 9 of Paper I) : 

w x v = TVS - VB. (73) 
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Using the property that the entropy and the Bernoulli 
constant are conserved along the stationary flow lines, 
and that the vorticity is zero at infinity, we deduce from 
Eq. ( |73| ) that the gradient of B is proportionnal to the 
entropy gradient and rewrite the vorticity as follows: 



7 



(74) 
(75) 



The increase of the temperature close to the accretor 
therefore implies that the contribution of the term V£> 
in Eq. ( f73| ) becomes negligible close to the accretor, like 
in BHL flows. 

The Bernoulli equation ( |36| ) is used to obtain an upper 
bound for the variation of the sound speed far from the 
accretor: 



(7 - 



< 



(7 

7- 



2GM 



(7 



(2GM\ 
\)M 2 V r J 



(M 2 - Ml) 



1 



,(76) 
(77) 



Using Eqs. ( pil| ) and ( |75| ) with optimal wavelength and 
neglecting compressibility effect, we obtain the following 
upper bound for the contribution of the region ahead of 
the accretor to the efficiency of the KHi : 



A 



KH 



< a- 



7 



1 



r A VS 



2+(7-l)M 2 oc 
< 0.07, 



log 



(78) 
(79) 



where r^/fA = 8, M c 



Thus both A^ T and 



0.6, and 7 = 5/3. 

are smaller than the esti- 



mated values of .Art (Eq. [40|) and Aku (Eq. |60|) in the 
vicinity of the accretor, although they might not be fully 
negligible. 

We conclude that these simulations , which might pos- 
sibly at first sight seem artificial, are an encouraging in- 
dication of the instabilities being generated by entropy 
gradients in the flow. They show that an instability exists 
which does not rely on the deformations of the shape of the 
shock surface, nor on the reflection of waves against the 
shock surface. Moreover, they suggest that the instability 
of the BHL flow is not an artifact due to the numerical 
treatment of the shock. 

These simulations should be followed by other simu- 
lations in order to understand better this instability. A 
first step would be to explore numerically the effect of the 
amplitude of the entropy gradient. The absence of shocks 
makes this flow easier to study by a global perturbation 
analysis in order to obtain conclusive statements about 
the precise onset of linear stability of the flow. Particular 
analytical solutions of such flows would be very useful in 
this respect. 



8. Conclusions 

Despite our lack of knowledge concerning both the shape 
of the shock surface, and the dependence of the shock dis- 
tance on the adiabatic index and Mach number, we are 
able to make a quantitative estimate of the entropy gra- 
dients produced by the shock. A priori, entropy gradients 
and vorticity are sources of linear instability through the 
Rayleigh-Taylor and the Kelvin-Helmholtz mechanisms. 
Their efficiency is estimated using a WKB like analysis, by 
integrating their local linear growth rate along a flow line 
between the shock and the accretor. The WKB criterion is 
only marginally satisfied, since the growth time is at best 
comparable to the advection time onto the accretor. This 
may cast doubts on the significance of our quantitative 
estimates, although the physical picture seems robust. If 
correct, this suggests that only large enough initial pertur- 
bations may reach non linear amplitudes when amplified 
by these mechanisms. 

It is striking that several features of the instability ob- 
served in numerical simulations (cf. the animations cur- 
rently available at http://www.mpa-garching.mpg.de/ 
-mor/bhla. html) would be explained naturally by these 
mechanisms: 

(i) the instability requires a shock: both the RTi and 
the KHi require the presence of a shock, which produces 
entropy gradients and vorticity. 

(ii) the instability is stronger when the shock is de- 
tached: the KHi adds to the destabilization of the flow if 
the shock is detached. 

(iii) the instability is stronger for high Mach numbers: 
the entropy gradient and vorticity produced by the shock 
increase with the Mach number. 

(iv) the instability is stronger for small accretors: the 
smaller the accretor, the longer the advection time, and 
the stronger the entropy gradients. 

(v) the instability is nonaxisymmetric: as is the RTi . 

(vi) the instability starts in the region of intermediate 
azimuthal angle (9 ~ 7r /2j, close to the accretor. this 
region coincides with the region of maximum entropy 
gradient and vorticity, if the shock is detached. 

If these mechanisms were responsible for the instabil- 
ity, our calculations further indicate that the efficiency of 
the linear instability becomes 

(i) independent of the Mach number if the kinetic en- 
ergy dominates the internal energy: A^^3>2/(7 — 1), 

(ii) independent of the accretor size if it is much 
smaller than the shock distance: r+ <C r n . Nevertheless, 
the range of timescales associated to the instability 
depends naturally on the size of the accretor. 

A higher efficiency of the linear instabilities would be 
reached if a feedback loop could be obtained: this is the 
subject of ongoing research. 
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Our study should invite the groups performing nu- 
merical simulations to follow carefully the evolution of 
entropy gradients and vorticity, since an overestimate of 
these quantities could lead to an artificially strong insta- 
bility. 

The simulations presented in Sect. (5] suggest that an 
absorbing accretor moving with a subsonic velocity in a 
gas with a non uniform entropy gives rise to an unstable 
accretion flow, resembling the BHL instability. This con- 
figuration might be a fruitful approach to understand the 
BHL instability better. 
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Appendix A: The Kelvin Helmholtz and 

Rayleigh Taylor instabilities in a 
compressible fluid 

A.l. Stratified sheared atmosphere 

We consider a stratified atmosphere with a density profile 
po(z), entropy So(z), sound speed co(z) in a vertical constant 
gravity go satisfying Eq. ( po| ) . The atmosphere is sheared in the 
y-direction with the velocity profile Vo(z). Perturbations are 
decomposed on a Fourier basis in the x, y directions, with their 
associated wavenumbers k x ,k y . Their growth rate is denoted 
by a. The linearized equations for perturbations lead to the 
following differential system on the perturbed vertical velocity 
v z and pressure p: 



obtain the following differential equation for the pressure per- 
turbation: 



^p 

dz 2 



, 9 log (kyVo 

dz 



dp 



+ 



^kyVo — icr y 



p = 



(A.3) 



The Orr-Sommerfeld equation governing the incompressible 
case (e.g. Drazin & Reid 1981) is recovered by taking the limit 
Co — ► oo. We solved this equation by a relaxation method be- 
tween two rigid walls, and obtained Figs. |^ and []. Different 
entropy profiles were tested, particularly the ones keeping a 
uniform Bernoulli constant with a length scale comparable to 
the vorticity length scale (as in an entropy-induced vorticity), 
with no significant departure f rom the case of uniform entropy. 
Indeed, we see from equation (A.3) that a non uniform entropy 
does not introduce significant changes compared to the case of 
uniform entropy, apart from the necessary spatial dependence 
of the sound speed co(z). 

A.3. Combined actions of shear and buoyancy 

The Richardson number Ri measures the ratio of buoyancy 
force to inertia for an incompressible fluid (see for example 
Chandrasekhar, 1961, §103): 



Ri : 



-go 



/dlogpN (dV \ 
V dz )\dz) 



rdVo- 

\ dz 



(A.4) 



It is also the squared ratio of the time scales associated to the 
RTi and KHi . In the classical case of an exponential density 
and a hyperbolic tan- velocity profile, the Kelvin-Helmholtz in- 
stability is stabilized for Ri > 1/4. 

Appendix B: Table of symbols used 



In Table B.l we present a list of the most important symbols 
used in this paper together with the equation number of their 
first occurence. 



d(ip Vz) 
dz 



dz 



k y V ' 



7-1 dS 



kyVo — ia 7 dz 
kyVo — icr x ' 



^kyVo — lay 



kl-kl 



(kyVo-ia) —9°-q^ 



kyVo — irr 
(ipoVz) 



(A.l) 



kyVo — icr 



- 9 4p 



(A.2) 



Taking the limit k x — > oo, k y = in Eqs. (A.l)- (|A,2j ), we 
recover the Brunt-Vaisala, frequency defined in Eq. (|21|). 

A.2. Compressible KHi without stratification 

The effect of com pres si bility on the linear KHi can be studied 
by solving Eqs. (|A.l|)-( [A~2j ) for g = 0. The pressure Po(z) 
across the flow is taken constant for equilibrium, but the en- 
tropy and the temperature may vary. 

Vo(z) being the unperturbed flow speed in the y-direction, we 
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